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Abstract 

In this paper we study a simple model of a purely excitatory neural network that, by construction, operates 
at a critical point. This model allows us to consider various markers of criticality and illustrate how they should 
perform in a finite-size system. By calculating the exact distribution of avalanche sizes we are able to show that, 
over a limited range of avalanche sizes which we precisely identify, the distribution has scale free properties but is 
not a power law. This suggests that it would be inappropriate to dismiss a system as not being critical purely based 
on an inability to rigorously fit a power law distribution as has been recently advocated. In assessing whether a 
system, especially a finite-size one, is critical it is thus important to consider other possible markers. We illustrate 
one of these by showing the divergence of susceptibility as the critical point of the system is approached. Finally, 
we provide evidence that power laws may underlie other observables of the system, that may be more amenable 
to robust experimental assessment. 



Introduction 

A number of in vitro and in vivo studies [IHl] have shown neuronal avalanches - cascades of neuronal firing - 
that may exhibit power law statistics in the relationship between avalanche size and occurrence. This has been 
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used as prima facie evidence that the brain may be operating near, or at, criticality [S]- In tm'n. these results 
have generated considerable interest because a brain at or near criticality would have maximum dynamic 
range [SHH] enabling it to optimally react and adapt to the dynamics of the surrounding environment (SUH] 
whilst maintaining balanced neuronal activity [10HT2]. Neuropathological states (e.g., epileptic seizures) 
could then be conceptualised as a breakdown of, or deviation from, the critical state, see [13j . for example. 
Furthermore, these findings have led to the notion that the brain may self-organise to a critical state [U], i.e., 
its dynamics would be driven towards the critical regime by some intrinsic mechanism and not be dependent 
on external tuning. In support of this view, Levina and colleagues |15j showed analytically and numerically 
that activity-dependent depressive synapses could lead to parameter-independent criticality. 

The interpretation that neuronal activity is poised at a critical state appears to be mostly phenomeno- 
logical whereby an analogy has been developed between the propagation of spikes in a neuronal network 
and models of percolation dynamics |16j or branching processes |17[|18j . Remarkable qualitative similarities 
between the statistical properties of neuronal activity and the above models have given credence to this 
analogy, however, the question remains as to whether it is justified. Indeed, various key assumptions under- 
lying percolation dynamics and branching processes are typically violated in the neuroscience domain. For 
example, full sampling, which is required in order to assess self- organised criticality, is unattainable even in 
the most simple in vitro scenario and yet it has been shown that sub-sampling can have profound effects 
on the distribution of the resulting observations 119] . On a related note, and more generally, the formal 
definition of a critical system as one operating at, or near, a second order (continuous) phase transition is 
problematic since the concept of phase transition applies to systems with infinite degrees of freedom |20| . 
Many neuroscience authors address this by appealing to the concept of finite size scaling and many published 
reports implicitly assume that distributions are power law with truncation to account for the so-called finite 
size effect. Typically, such reports adopt an approach whereby (a) scale invarianee is assessed through finite 
size sealing analysis, confirming that upon resealing the event size, the curves collapse to a power law with 
truncation at system size (but see below regarding the definition of system size); (b) the parameters of 
statistical models are estimated, typically over a range of event size values that are rarely justified; and (c) 
the best model is determined by model selection, in which power law and exponentially truncated power 
law are compared to alternatives such as exponential, lognormal and gamma distributions, see |21| for a 
typical example. Whilst greater rigour in the statistical treatment of the assessment of the presence of power 
laws has been attained following Clauset and colleagues' infiuential paper [35], what seems to be lacking is 
a rigorous treatment as to why a power law should be assumed to begin with. Although this question is 
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particularly pertinent to the neurosciences, it should be noted that similar questions remain open in the field 
of percolation theory (e.g., namely: (i) how does the critical transition behaviour emerge from the 

behaviour of large finite systems and what are the features of the transition? (ii) what is the location of the 
scaling window in which to determine the critical parameters? 
This paper specifically seeks to address the following questions: 

1. Assuming that the whole brain, or indeed a region of interest defined by where data can be obtained, is 
operating at criticality, can we reasonably expect power law statistics in neural data coming from a very 
small (possibly sub-sampled) subset of the system? If not, what would be the expected distribution? 
Sornette j25j states that the Gamma distribution is "found in critical phenomena in the presence of 
a finite size effect or at a finite distance from the critical point" . Jensen |26| claims that finite-size 
systems often show an exponential cut-off below the system size. However, we are not aware of any 
study in which the distribution of event sizes in a finite-size system set to operate at a critical regime 
has been investigated. 

2. In a finite-size system, is it reasonable/possible to perform a robust statistical assessment of power law 
statistics? Even the application of a rigorous model selection approach will lead to different results 
depending on the choice of the range of event sizes and the number of samples being considered |27j . 
The issue of range selection is of particular interest. Whilst the notion of system size is clear in models 
of criticality such as the Abelian sandpile where (i) there is full sampling, (ii) the number of sites is 
finite, and (iii) there is dissipation at the edges, system size is much less obvious where re-entrant 
connections exist, making it possible, in principle, for avalanches to be of infinite size. Here, the 
counting measure which leads to the definition of an avalanche is important. Counting the number 
of neurones involved in an avalanche will lead to a clearly defined system size, whereas counting the 
total number of activations - the de facto standard, e.g., [TIIlfTSlI^ - will not. Furthermore, it should 
also be noted that the presence of re-entrant connections invalidates the standard theory of branching 
processes |18] . and makes a rigorous determination of the branching parameter a problematic if not 
impossible, e.g., in the presence of avalanches merging. 

3. Are there other markers of criticality that may be more amenable to characterisation and that should 
be considered instead of, or in addition to, the statistics of event sizes? The need for such markers 
in neuroscience has been recognised (see j27j for example) and a number of studies have investigated 
long-range temporal correlations (power-law decay of the autocorrelation function) in amplitude fluc- 
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tuations 129' and in intcr-burst intervals }30[I31| . However, a theoretical account of how those may 
relate to one another is lacking (although see the recent work in [32] )■ Other markers of criticality (or 
markers of transitions) have been associated with critical physical systems, e.g., divergence of suscep- 
tibility and slowing of the recovery from perturbations near the critical point |25j . however, we are not 
aware of any theoretical or empirical study investigating them in a neuroscience context. 

One way to address these questions more rigorously is to use simplified but therefore more tractable 
conceptual models (e.g., j33|). In this paper, we use a model of a purely excitatory neuronal system 
with simple stochastic neuronal dynamics which can be tuned to operate at, or near, a second order phase 
transition (specifically, a transcritical bifurcation). The simplicity of the model allows us to analytically 
calculate the exact distribution of avalanche sizes, which we confirm through simulations of the system's 
dynamics. We study our model at the critical point and compare our exact distribution to the explicit but 
approximate solution proposed by Kessler |34j in an analogous problem of modelling disease dynamics. We 
confirm that Kessler's approximate solution converges to our exact result. Importantly, we show that, in the 
proposed finite-size system, this distribution is not a power law, thus highlighting the necessity of considering 
other markers of criticality. We thus analyse two potential markers of criticality: (i) the divergence of 
susceptibility that arises in the model as we approach the critical point, (ii) the slowing down of the recovery 
time from small disturbances as the system approaches the critical point. Finally, we speculate on a sufficient 
but not necessary condition under which our exact distribution could converge to a true power law in the 
limit of the system size. 

The stochastic model 

We start from the stochastic model of Benayoun et al. [TU] which we simplify to the most trivial of models. 
A fully connected network of N neurones is considered with purely excitatory connections (as opposed to the 
excitatory and inhibitory networks considered in [lOj). Within the network, neurones are considered to be 
either quiescent (Q) or active (A). Taking a small time step dt and allow rfi — > the transition probabilities 
between the two states are then given by: 

P{Q ^ A, in time dt) ^ f {si(t)) dt 
P {A-^ Q, in time dt) ~ adt 

where Si(t) = -ff-aj{t) + hi is the input to the neurone. Here / is an activation function, hi is an optional 
external input, Wij is the connection strength from neurone i to neurone j and aj{t) = 1 if neurone j is 
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active at time t and zero otherwise, a is the de-aetivation rate and therefore controls the refractory period 
of the neurone. 

To aUow tractabihty, we further make the following simplifications: 

1. We assume that all synaptic weightings are equal (wij = w). 

2. We assume there is no external input. The driven case will be explored theoretically and empirically 
in a companion manuscript. 

3. We assume the linear identity activation function f{x) = x. Although it is more common to use 
sigmoid activation functions we note that the identity function can just be thought of as a suitably 
scaled tanh function over the desired range. 

As the network is fully connected we can write the following mean field equation for active neurones: 

dA wA ^ , wA 

— = Q -aA= [N - A)- aA, 

where we have appealed to the fact that the system is closed and thus A + Q ^ N. This ODE is analogous to 
the much studied |35j susceptible — infectious — > susceptible (SIS) model used in mathematical epidemiology 
and we can appeal to some of the known results in studying its behaviour. Primarily we can use simple 
stability analysis. The non-zero steady state is given by A* = N{1 — a/w). Setting g{A) = this 
equilibrium point is stable if g'{A*) < 0. Thus 

g'iA) ^{w^a)- 2w^ ^ g'{A*) = {w - a) - 2w ^^^ "J^^"^^ ^a-w. 

Borrowing from epidemiology we define the threshold Rq = ^- li Rq > I we see that g'{A*) = a — w < 
and the non-zero steady state is stable. Figure [T] illustrates the differing behaviour of the solution to the 
above ODE for Rq < 1 (sub-critical), Rq = 1 (critical) and i?o > 1 (super-critical). 

Firing neurones and avalanches 

Instead of focussing on the average activity level across the network we will instead look at the size distri- 
bution of firing neurones following one firing event. To do this we begin with a fully quiescent network and 
initially activate just one neurone. We then record the total number of neurones that fire (the number of 
quiescent to active transitions) until the network returns to the fully quiescent state. We use this process 
of sequential activation as our definition of an avalanche and our main interest is the distribution of the 
avalanche sizes. Unfortunately, the simple ODE approach will not provide us with this distribution. To 
calculate this distribution, we use the semi-analytic approach described in the following section. 
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Figure 1: Activity in the different regimes. Plot of the solution to the ODE for N = 800 and three 
different regimes where Rq was set to 0.5 (blue), 1.0 (green) and 2.0 (red). Initially we activated 25% of the 
network. 

Tree approach to the avalanche distribution 

We begin by defining qi as the probability the next transition is a recovery (from A to Q) given i active 
neurones (i > 0). The probability the next transition is an activation is then 1 ^ qi and we have: 

_ aN _ N 

* ^ w{N -i) + aN ^ R(i{N -i) + N' 
w{N-i) _ Ro{N-i) 
^ ^ w{N -i) + aN ^ Ro{N -i) + N' 

In order to calculate the avalanche size distribution we adopt a recursive approach. We begin by considering 
the process unfolding in a tree like manner with 1 initially active neurone. The tree can be divided into levels 
based on how the process is unfolding. Each level is itself subdivided into two containing the possible number 
of active neurones before and after a transition has occurred. Level k contains firstly the number of active 
neurones after 2k transitions and secondly all possible active neurones after 2fc-|- 1 transitions. The recursive 
tree approach relates the probability of transition between and within these levels to the final avalanche size. 
Figure [2] illustrates the first 3 levels of this process. To continue we define p^, as the probability of having i 
active neurones on level k with i = 0,1,2, N and k G Nq. On the first level {k = 0) we immediately see 
that Pq = 1, Pq = 1 — qi and p^ = qi. To proceed we will consider the probability of having an odd number 
of active neurones on an arbitrary level. First we note the following within level recurrence 

Pk = Pk"'^ (1 - 92j-i) + pl^^^q-ij+i [j > 1). 
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We can now use these in our calculation of pf^^i and obtain the following between level recurrence 

2j+l 2j /-, N , 2j+2 

=Pk~^ (1 - 92j-i) (1 - q2j) + pI^^^ q2j+i (1 - q2])+pi^^^ (1 - 92j+l) q2j+2 + pl^^^q2j+3q2j+2 

=Pk^^ (1 - fej-i) (1 - 92j-i) (92J+1 (1 - q2j) + (1 - q2]+i) q2]+2) + Pk'^^q2j+3q2j+2- 

This newly derived recurrence relation provides a consistent set of equations between levels in terms 
of odd numbers of active neurones. This algebraic manipulation was necessary to obtain a self-consistent 
system of equations. Letting rj(k) = p^^^^ we have: 

rj(fc + 1) ^rj_i{k) (1 - g2j-i) (1 - q2j) + rj{k) {q2j+i (1 - q2j) + (1 - 92j+i) 92j+2) + fj+i (fc)92j+392j+2, 

where we must make suitable modifications for systems with only an individual neurone active and also near 
and at full activation. For a single active neurone the equation is simply given as 

ro(fc + 1) =ro(fc)(l - qi)q2 + ri{k)q3q2. 
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For full (or near full) activation the equations depend on the system size, A^. Defining N = {N — 2)/2 if N 
even and N = {N — l)/2 for N odd we then obtain 

+ I' w-i(^) (1 - <l2N-i) (1 - 12n) + ^n(^) {I2N+1 (1 - 12n) + (1 - I2N+1) 92JV+2) ' if N even, 
^ \'"7v-i(^)(l-92W-i)(l-92Jv) + ''w(^)92W+i(l-92Jv). otherwise. 



We now define: 



/ roik) \ 



r{k) = 



We can now write r(fc + 1) 



matrix: 



/ (l-(7i)(72 9392 



A • r(k) where, for N even, matrix A is given by the following tridiagonal 



bj '^j 



\ 



(1 - <l2N-l) (1 - 92JV) <l2N+l (1 - 12n) + ^2N +2 (l " I2N +l) ' 



where 



(1 - q2j-i) (1 - q2j) 



Cj = q2j+i (1 - q2]) + q2j+2 (1 - 92^+1 ) 

dj = (j2j+3'72i+2 and j = 1 : {N - 1). 

On the /c*^ level of the tree, the probability of only 1 neurone being active is given by = rQ(fc). We can 
then calculate the probability of zero active neurones after k firings as qiro{k)] this is thus the probability, 
P{k + 1), of having avalanches of size k + 1 since initially one neurone is active. To calculate the distribution 
we implemented the recursive method of calculation in the MATLAB® environment. 



Simulations of neuronal avalanches 

In order to check the validity of our method, we performed simulations of the firing neurones using the 
Gillespie algorithm |36j . Due to the network being fully connected the algorithm is relatively straightforward: 



• As earlier let A be the number of active neurones in the network {Q the number of quiescent). Given 
that an individual neurone becomes quiescent at rate a then the total rate of (Active — >■ Quiescent) 
transitions is given by — Aa. Similarly the total rate of (Quiescent — > Active) transitions is given 
by r,a = fis^)Q^f{s,){N-A). 

• Let r = Taq + Tqa and generate a timestep dt from an exponential distribution of rate r. 

• Generate a random number n between and 1. If n < ^ an active neurone turns quiescent, otherwise 
a quiescent neurone is activated (fires) . This event is said to occur at time t + dt and the network is 
updated accordingly. 

Exact solution compared to simulation 

Values of the threshold, Rq, were chosen less than, equal to and finally above 1. We will refer to these regimes 
as subcritical, critical and supercritical respectively. Figure [3] illustrates the, as expected, good agreement 
between the simulations and the exact result for the three different regimes of Rq. 
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Figure 3: Avalanche distributions. Results from the simulations of the avalanche distributions for the 
subcritical {Rq < 1, blue), critical {Rq = 1, green) and supercritical {Rq > 1, red) regimes for a network of 
size N = 800. For each regime 2, 000, 000, 000 avalanches were simulated. The corresponding exact solutions 
are shown in black. 

A closed solution 

In [34], Kessler proposed a closed solution to the analogous susceptible- infected-susceptible (SIS) problem 
where he was interested in the number of infections (including reinfections) occurring over the course of an 



9 



epidemic. For small avalanche sizes where the number of infectives is negligible compared to the network 
size the transition probabilities can be approximated as 

N 1 

1 - <?z 



Ro{N-i) + N i?o + l 
RoiN-t) Ro 



Ro{N-i) + N Ro + l' 

In the critical regime Rq = 1, the problem reduces to calculating the distribution of first passage times of 
a random walk with equal transition probabilities. Thus for avalanche sizes in the range, 1 <C n <C x/iV, 
Kessler [33] gave the following distribution based on Stirling's approximation 



P{n) 



2n - 2\ /2n - 2 
n — 1 7 \ n 



1) 



We note however that the range over which the distribution can be shown to be a power law is rather limited 
and for small networks will not hold. Using the theory of random walks and a Fokker-Planck approximation 
Kessler also derived the following closed solution to the probability distribution of infections in the critical 
regime (i?o = 1) for larger sizes: 

1 3 

P{n) = , exp(n/2iV) sinh" ^ (n/N) {n > 1) (2) 

Figure 2] plots this approximation against our exact solution for a network of size N = 800. To more 
formally assess the convergence of the approximate solution to that of our exact solution we considered the 
probabilities of avalanches from size N/IQ to 20N and measured the difference between the distributions 
using two different metrics. Letting Pe{n) be the exact probability of an avalanche of size n and Pk{n) be 
the Kessler approximation to this we first considered the standard mean-square error given by 

^ 20N 

Error{N) ^ - ^ {Pe{n) - Pk{n)f where R = 2QN - N/IO + 1. 

n=N/10 

Secondly we considered a more stringent measure of the error by looking at the supremum of difference 
between the same range of avalanches 

Err or {N) = sup \Pe{n) - Pk{n)\. 

n 

Figure [5] illustrates the two errors for increasing network size and both show how the proposed closed solution 
is indeed converging to that of the exact. 
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Figure 4: Closed solution versus exact. Plot of the closed solution (red solid line) versus the exact 
solution (black dashed line) for a network of size N = 800 operating in the critical regime. 




Network size {N) Network size (N) 

(a) Mean square error (b) Suprcmum error 

Figure 5: Convergence of closed solution to exact, (a) Here the mean square error is given by the blue 
line and 0(-/V^) and 0{N^) convergence represented by the black and red lines respectively, (b) Here the 
suprcmum error is given by the blue line and 0(A^^) and 0(A^^) convergence represented by the black and 
red lines respectively. 

Scale-free behaviour in the Rq = 1 regime 

Whilst Equation [1] gives a power law, this equation only holds over a limited range. Equation [21 in turn, is 
neither a power law nor a truncated power law. Here, we assess the range over which the distribution of sizes 
can be said to exhibit scale-free behaviour. For a rigorous assessment of this range, we employ a subset of the 
model selection approach described by Clauset and colleagues [22j. Specifically, we consider 100,000 of the 
simulated avalanches described earlier and fit a truncated power law distribution of the form P{x) = Cx~" 
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to avalanches up to size Xmax = {the choice of this upper bound will be justified in the following 

section) by using the maximum likelihood method (here C is a normalising constant to keep the sum of the 
distribution between [xmim Xmax] equal to 1). We do this by finding values of a and Xmin that maximise the 
probability of obtaining our simulated avalanches given the fitted distribution. Next we randomly generate 
1000 data sets from the fitted distribution and compute the difference between these synthetic data sets and 
the fitted form (using the Kolmogorov-Smirnov statistic). Similarly we compute the difference between our 
simulations and the fitted power law. The p- value is then calculated as the proportion of synthetic data sets 
that arc further away from the theoretical distributions than our simulations. As per [22], the hypothesis 
(that the data comes from a power law) is rejected if the p-value is less than 0.1. Note that in the model 
selection approach, should the hypothesis not be rejected, then one should test alternative models and use 
an information criterion to identify the best model. However, our focus here is purely on assessing whether 
our distribution can be said to behave like a power law distribution (we know it is not actually a power 
law) and therefore alternative models were not tested. With 100, 000 avalanches we obtained a p-value of 
0.382 leading us not to reject the hypothesis that the distribution was power law (see figure [6]). Since the 
distribution is not a power law, we would expect that upon considering a larger number of avalanches, this 
hypothesis should be rejected [21]. Indeed, using data from 1, 000, 000 avalanches yielded a p-value of 0, i.e., 
the truncated power law is not an appropriate model for the distribution. 



Avalanche size (71) 

(a) PDF 




Avalanche size [n) 

(b) CDF 



Figure 6: Fitted distributions. Out of 100, 000 of the observed avalanches we fit the 98, 833 whose size 
was less than j^N . (a) The fitted probability distribution function (black line) fitted over the simulated 
avalanche distribution (blue dots), (b) The fitted cumulative distribution function (black line) fitted over 
the simulated avalanche distribution (red dots). 



The fact that the truncated power law was a plausible fit for the lesser number of avalanches (note that 
100, 000 is of the same order of magnitude as the number of avalanches typically reported in in vitro or in vivo 
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studies of neuronal avalanches) is indicative of the partial scale-free behaviour the model exhibits. The appeal 
of the concept of critical brain is that the critical regime is the one in which long-range correlations keep the 
system poised between too highly correlated states of no behavioural value and too weakly correlated states 
that prevent information flow |37j . Thus, the actual nature of the distribution of the avalanche size matters 
less than any indication of the presence of long range correlations. In other words, neuronal avalanches need 
not precisely follow a power law, they just need to exhibit similar behaviour. It is important to appreciate 
this distinction. As the exact solution to the distribution of avalanche sizes is known, we can then compare it 
visually with a fit of a truncated power law over avalanche sizes from jqN to jqN . This is done in Figure [7] 
which confirms that over a limited range of sizes the distribution is well approximated by a power law. 




10" I ' — ' — ^ ' — ' — ' — 

10° lo' 10" lo' 10* lo' 

Avalanche size (n) 

Figure 7: Power law fit of the exact solution. Main: plot of the truncated power law (red) fitted 
over the entire range of the exact distribution (black). Inset: Fitted power law and exact distribution in the 
range {^N,^N]. 

Origin of the distribution's truncation 

The fact that we have an exact form for the distribution allows us to make further important observations 
about some of its characteristics. Here, we explore the origin of the distribution's truncation. Let Aq, Ai, 

, be the eigenvalues of A with the corresponding eigenvectors mq, ui, , Ujy. The initial 

condition can then be given as r(0) = cqMo + ciui + + CfjU^. Using the property A.Uj — XjUj we then 

obtain r(fc) = cqAqUo + ciAf ui + + c^A^Ujy- This calculation leads to the probability of an avalanche 
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being of size n being: 

N 

P{n)^q,Y.d.X:, (3) 

1=1 

where qi is the probabihty that the next transition is a recovery (from A to Q) given 1 active neurone 
(as defined carher), A^s are the eigenvalues of the transition matrix A and djS are specified by the 
eigenvectors of the transition matrix and the initial conditions. We note that the earlier cciuation, 

r(0) = cquq + ciui + + CffUfj^ can be solved to obtain a. Using this, we can then calculate di as 

the first entry of the vector CiUi. Equation [31 which is exact, thus demonstrates that the distribution of 
avalanche sizes is a linear combination of exponentials. 



Assuming the lead eigenvalue is denoted by Ai, then for all i, Ai < Ai and we have 



N 



P{n) _ ^ diX. 



di \\i di [Xj ^■•■^ di I Ai 



Taking the limit as n increases we find 

n-i-oo qidiXi 

Hence P{n) ^ qidiX^ and for larger avalanche sizes we have the leading eigenvalue dominating thus giving 
the exponential cutoff observed. We illustrate this convergence in Figure [S] where we plot the exact avalanche 
distribution, P{n), against qirfiA". This figure also illustrates that the leading eigenvalue begins to dominate 
for avalanches just over the system size. It is for this reason that we chose an upper bound of ^ when 
fitting a power law to the distribution of avalanche sizes in the previous section. 

Power law distribution for large systems 

Whilst we are unable to show analytically that in the limit of the system size the distribution converges to 
a power law, we can conjecture that this is the case. We motivate this by considering a hypothetical system 
that can again be characterised by its transition matrix and Equation [3] for the probability distribution. 
Such a distribution could converge to a true power law under two important conditions: 

1. the eigenvalues of the transition matrix A arc well approximated by a geometric distribution, i.e. they 
are in the form Ai = A'e^^', 
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Figure 8: Exponential cutoff. Exact avalanche distribution (black line), plotted against a distribution 
assuming only the leading eigenvalue is non-zero (blue line). Avalanches greater than the system size, 
N = 800, appear after the dashed line. 

2. the constants di in Equation [3] are well approximated by di = Li'^, 

where K, fi, L and q can be inferred via a numerical fitting procedure. In such a scenario. Equation [3] can 
be rewritten to give 

TV 

P{n) = CY,^%en-\ (4) 

4=1 

where C is a given constant. In the limit of an infinite network size we then have 

oo 

P{n)^CY,^'^{en-'■ (5) 

i=l 

While P{ti) can be found based on standard mathematical arguments, we have chosen to use results derived 
in the context of the Z-transform. The standard results for integer values of q give 

oo 

Y^^''^-^ = (-l)''D'^{^), (6) 
^-^ z — 1 

i=l 

where D is an operator such that D{f{z)) = z^^^^^. For a fixed integer value of q, an approximation for 
P{n) can be obtained by simply applying the operator as many times as necessary and then substituting 
z = e'^". For q ~ 1 for example, P{n) cx ^^f,',_-^yi which for small values of fi is well approximated hy -i^^. 

These results only hold for integer values of q so an alternative approach is to approximate the sum for 
P{n) in terms of an integral. Taking into account the special form for the eigenvalues and constants, P{n) 
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can be approximated as follows: 

P(n) =C^z'?(e^")-'~C / x'^e-'^^'^dx. (7) 

The latter integral can be interpreted as a Laplace transform of and thus yields 

PW^C^fcil^. (8) 

It is worth noting that this result is consistent with that obtained for integer values of q. 

For a simple empirical verification of this conjecture, we determined the values of K, fi, L and q in the 
above conditions that fitted the first and twentieth eigenvalues and d constants of the exact distribution for 
a network of size = 800 (sec figures 9(a) and |9(b)] ) and compared the resulting probability distribution 
with the exact distribution. As shown by figure |9(c)[ there is remarkable agreement between both curves 
over a broad range of values, including the range [jqN, jqN] over which a power law like behaviour was 
established earlier (see figure [T]). This result clearly illustrates the dominance of the larger eigenvalues and, 
given that the fitted distribution converges to a power law, gives support to the conjecture that the exact 
distribution would do so in the limit of an infinite network. 



Other markers of criticality 

Since the distribution of avalanche sizes in the finite-size critical system does not necessarily follow a true 
power law, the application of robust statistical testing in experimental conditions could well lead to rejecting 
the hypothesis that the data may come from a system operating in the critical regime. Therefore, in this 
section, we consider two experimentally testable markers of criticality: critical slowing down and divergence 
of susceptibility. We will define those concepts below but first we briefly summarise Van Kampen's system 
size expansion [38] which we use to illustrate those markers on our system. 

System size expansion 

For generality we now assume that each neurone receives a constant external input and that the activation 
function can take forms other than the simple identity function. We define the probability that the number 
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(a) distribution of eigenvalues 



(b) distribution of constants di 



(c) distribution of avalanche sizes 



Figure 9: Possible origin of the power law for large systems, (a) Actual distribution of eigenvalues 
Xi (black crosses) along with fitted distribution (blue line), (b) Actual distribution of constants di (black 
crosses) along with fitted distribution (blue line), (c) Exact distribution of avalanche sizes (black) along 
with distribution resulting from fitted distributions of and di (blue). All plots are for a network of size 
N ~ 800 operating at criticality. 



of neurones active at time f is n as P„(i). Then the master equation can be given as 

dPn (t) 



dt 



=a (n + 1) P„+i {t) - 
anPn {t) + 



f{^ + h) {N-n)P^{t). 

The idea is to now model the number of active neurones as the sum of a deterministic component scaled by 
N and a stochastic perturbation scaled by \fN , i.e., 



nit) = Nn{t)+N^^{t). 
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Full details of the system size expansion are given in the Appendix, but what we obtain is the following set 
of equations for /i (which is the solution to the mean field equation of the proportion of active neurones), 
(5) (the expected value of the fluctuations) and ct^ — ^^^) — (^)^ (the variance of the fluctuations) 

§^ = -«M+(1-m)/, (9) 



d{0 
dt 

dt 



a 



+ f-wf'il-^i)){0, (10) 



-2 (a + f- wf'il - m)) (^2) + (afi + (1 - /i)/) . (11) 



Here / = fiwfi + h) and /' = f'{w^ + h). These equations, in turn, give the following equations for the 
mean, A, and variance, A„, of the number of active neurones 

A = N 1.1 + iV^2 — (assuming we know the initial number of active neurones), (12) 

A,^N{a^). (13) 

Critical slowing down 

In dynamical systems theory, a number of bifurcations, including the transcritical bifurcation in our system, 
involve the dominant eigenvalue characterising the rates of changes around the equilibrium crossing zero. 
As a consequence, the characteristic return time to the equilibrium following a perturbation increases when 
the threshold is approached [39]. This increases has led to the notion of critical slowing down as a marker 
of critical transitions |40| . Here, we illustrate the critical slowing down of our model with the analytic 
derivation of the rate of convergence to the steady state (this derivation has been previously shown by [41]). 
We first begin by calculating the analytic solution to Equation [9] for the percentage of active neurones. We 
again consider the case where / is the identity function and can thus write 

^ = -a^l + (1 - ^l)f{w^l + h) = -a^i + (i - ^l){w^l + h). (14) 

Assuming zero external input (/i = 0), we have 

— = —afi + (1 — fJ.)(wii + h) — —afj, + (1 — ^)wfi. (15) 

We are interested in the solution of this equation and consider the result for different values of a. Firstly we 
consider a ^ w. In this case we have 

d Ijj 

— = — a/i + (1 — fj,)w^ = /i(w — w/i — a). (16) 
ot 
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Integrating this using separation of variables and the initial condition fi(0) — ^iq, we find 

M = -n u where A = . 17 

The solution to this depends on whether a < w ot a > w (Rq > 1 and i?o < 1 respectively). If a < w then 
as t — > oo, /I — > li a > w then as i oo, /i 0. Note that in both cases, convergence of the number 

of active neurones to the steady state solution is exponential. 

Now we consider the solution when a = i.e., the critical regime. 



— = —afj, + (1 — /i)a/i — ~afi fi(t) 



Thus as i — !■ oo we find — > 0. However, unlike for Rq ^ 1, convergence to the steady state exhibits a 
power law dependence on time [41) . 

Divergence of susceptibility 

A correlate of the phenomenon of critical slowing down is that of the divergence of susceptibility of the 
system as the system approaches the bifurcation |40] . In this section, we investigate the behaviour of the 
equation for the variance. For simplicity, we consider again the case of the identity activation function and 
a non-driven system h = 0. First we use Equation [TT] to calculate the variance in the percentage of active 
neurones: 

o 2 

-^ = -2(a + f- wf (1 - /i)) + (afi + {1 - ^i) f 



dt 

= —2 (a + wii + /i — (1 — /i)) (T^ + (a/i + {1 — jj) (w/i + h)) 
= — 2 (a + wfi — {1 — /i)) + {afj. + (1 — /i) wfi) 



Setting this equal to zero and rearranging we obtain 



(a/i + (1 - /x) w^i) _ (/i + (1 - /z) Ron) 



2{a + wfi-w'^ {1- /i)) 2(1 + Rofi - Rqw (1 - fi)) ' 

Here we note that unlike the equation for fi where there was only the single bifurcation parameter Rq, we 
now have the additional dependence on w. To maintain consistency with earlier results, we now set w = 1 
to obtain 



a if a < 1 (i?o > 1), 
lim (7^(0 =<(! ifa = l (i?o = l), 

otherwise (i?o < !)• 
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Using Equation [T3] we obtain 



N 



if Ro > 1, 



t— >-oo 



if i?o = 1, 
otherwise {Rq < 1). 



Figure [TOl illustrates the jump to a non-zero steady state when the critical value i?o = 1 is approached from 
below, and the divergence in variance when it is approached from above. 



800 




Figure 10: Divergence of susceptibility. Analytic result for the steady state of the variance as Rq 
approaches 1 in a network of size = 800. Results only provided down to a = 2/3 for clarity. 

Here it should be noted that any finite-size network has a zero absorbing state so that eventually all 
activity will die out irrespective of the value of Rq. However, it has been shown that the ODE limit is 
a valid approximation to the solution of the master equation for reasonably sized systems with values of 
Rq greater than 1 and only over a finite time horizon (see j42j for further discussion). Defining the true 
(i.e., calculated directly from the master equation for P(n)) expected value of active neurones at time t as 
A{t), the convergence of the ODE approximation for A{t) given by Equation [T^] is such that for any t > 
hm \A{t)~A{t)\=Om- 

N^OG 



Discussion 

Over the last decade or so, the search for evidence that the brain may be a critical system has been the focus 
of much research. This is because it is thought that a critical brain would benefit from maximised dynamic 
range of processing, fidelity of information transmission and information capacity |44) . Whilst support for 
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the critical brain hypothesis has emerged from comparing brain dynamics at various scales with the dynamics 
of physical systems at criticality (e.g., [inHHHH]), in this paper, we focus on the important body of work 
that has relied on characterising power laws in the distributions of size of neuronal avalanches [71l28j . Our 
focus on this scale is motivated by empirical considerations regarding how one can go about demonstrating 
the above functional properties. Shew and Plenz |44| remark that any research strategy to test whether 
these properties are optimal near criticality will have to achieve two criteria: a means of altering the overall 
balance of interactions between neurones and a means of assessing how close to criticality the cortex is 
operating. As argued by these authors, the study of neuronal avalanches offers the greatest likelihood of 
achieving those two criteria. 

The importance of a robust assessment of the statistical properties of the avalanche size is therefore 
two-fold; on the one hand, it is about ascertaining the extent to which the system being studied has the 
statistical properties expected of a system operating at, or near, criticality; on the other hand, it is about 
being able to confirm that a manipulation/perturbation of the system aimed to push the system away 
from this critical regime has been effective. This consideration therefore puts a lot of importance on the 
description of the statistics one should expect in such a system. In the current literature, the assumption 
of the distribution of avalanche sizes taking a power law functional form relics on an analogy between the 
propagation of spikes in a neuronal network and models of percolation dynamics or branching processes for 
which exact power laws have been demonstrated in the limit of system size. As a result of the importance of 
having a robust assessment of the expected presence of a power law, greater emphasis has recently been put 
on using a sound statistical testing framework, e.g., |22| . Whilst we are unaware of any study in which the 
criticality hypothesis was rejected due to failure of rigorous statistical testing (which we suspect is due to the 
necessarily small number of observations, as we will argue below), there is clear evidence that many authors 
are now using the methods of Clauset et al. [22] to confirm the criticality of their experimental findings, 
e.g., [iniHIlUT]. As a result, we feel that it is all the more important to confirm that the assumed power 
law functional form is indeed a sensible representation of what one should expect in in vivo and in vitro 
recordings, which, unlike the physical systems considered when deriving the power law statistics, are finite- 
size systems. The aim of the paper was therefore to consider a model of neuronal dynamics that would be 
simple enough to allow the derivation of analytical or semi-analytical results whilst (i) giving us a handle on 
the parameter controlling the fundamental principle thought to underlie criticality in the brain, namely, the 
balancing between processes that enhance and suppress activity (note that we are intentionally not referring 
to excitation and/or inhibition - we will return to this below) and (ii) allowing us to determine its distribution 
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of avalanche sizes when operating in the critical regime. Note that because we are using a finite-size system, 
we are appealing to a normal form of standard bifurcation, here, a transcritical bifurcation, because it 
embodies all that needs to be known about the 'critical' transition (Sornettc, private communication). 

Our semi-analytic derivation of the true distribution of avalanche sizes in a finite-size system suggests 
that, even though it is approximately scale free over a limited range, the distribution is not a true power law. 
First, this has important implications for the interpretation of results from a robust statistical assessment of 
the distribution. Indeed, as has been discussed by Klaus and Plenz [2 1 ' , with a large number of samples, any 
distribution that deviates from the expected distribution by more than noise due to sampling, will eventually 
yield a p-value such that the power law hypothesis will be rejected, thus leading to the potentially incorrect 
conclusion that the system is not critical. This is the case in our scenario where using 10^ avalanches 
lead to a rejection of the criticality hypothesis even though the system is tuned to the critical regime. In 
contrast, with 10^ avalanches (which is consistent with empirical observations), a p-value above threshold 
leads to not rejecting the hypothesis that the distribution is a power law even though we established it is 
not one^. This finding therefore provides an important counterpart to the analytical results of Touboul and 
colleagues [27] who showed that thresholded stochastic processes could generically yield apparent power laws 
that only stringent statistical testing will reject. Whilst the stringent testing will reject the hypothesis of 
criticality for a system that is not necessarily critical, it may also reject the hypothesis of criticality for a 
system that is critical only because the actual distribution is not actually a power law. This ambiguity of 
the avalanche distribution in the finite-size system therefore requires that one should carefully consider to 
what fundamental property the idea of a critical brain actually appeals to. We suggest that the key appeal 
is that the brain can exhibit long-range correlations between neurones without it ever experiencing an over 
saturation of activity or long periods of inactivity. It then follows that the importance is not in the exact 
distribution obtained but in the approximately scale-free behaviour it exhibits. In turn, this highlights the 
importance of looking at other markers of criticality (which we will discuss below). 

Another important result of this work is to provide the beginning of a mechanistic explanation for an often 
alluded to (e.g., [SO]) but never properly treated (as far as we are aware) observation that whereas avalanches 
in a critical system with re-entrant connections could in principle be arbitrarily long, and certainly, exceeding 
the number of recording sites, neuronal avalanches in in vitro or in vivo systems (and many computational 
models of self-organised criticality) often show a cut-off at the number of sites. Our work suggests that the 
lead eigenvalue of the transition matrix between states fully determine the location of this cut-off which 
turns out indeed to be at about the system size, even if avalanches of up to 20 times the system size can be 
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observed. This finding therefore provides some justification for setting, or accepting, a bound within which 
to apply a Clauset-type methodology (we note that various reports use different ranges, e.g., 80% of system 
size in |15] . roughly system size in |50j). It is worth remembering that the number of recording sites can 
have profound implications on the nature of the distribution observed [19j . 

In addition to providing results on the distribution of avalanche sizes, we also sought to explore other 
potential markers of criticality. We provided results on two other markers of criticality - critical slowing 
down and divergence of susceptibility - both of which again follow from a dynamical systems appreciation 
of a critical bifurcation, i.e., the behaviour of a system whose lead eigenvalue crosses zero. The appeal of 
those markers, which have been documented in many other natural processes, e.g., |40|j51j . but seldom at the 
mesoscopic brain leveP, see |52| for a rare example, is that (a) they strengthen the assessment of the system 
being critical and (b) may contribute to achieving the second criterion of Shew and Plenz [44]. Although 
the authors are not in a position to provide explicit recommendations for an experimental design, we believe 
that these markers are amenable to robust experimentation, e.g., through pharmacological manipulation. 

Whilst we hope we have convinced the reader of the potential importance of these findings, we also need 
to recognise that the very simplicity that makes analytical work possible does also raise questions regarding 
how physiologically plausible such a model is and therefore whether its conclusions should be expected to 
hold. Below, we address a few of the points worthy of further consideration. 

Validity of a purely excitatory network 

In this paper, we have used a purely excitatory neuronal model. This not only simplifies the system but is 
also an important characteristic of the brain during early development. Experimental results have shown that 
during early development, before birth, GABAergic neurones (i.e. neurones which will later be inhibitory) 
have a depolarising effect on their post-synaptic neighbours [55H55] . Thus, our model might be considered 
as representative of early development. Power law statistics have been observed in early development at 
a time when networks are thought to be purely excitatory |30[l56j . It should be noted that this approach 
has the benefit of casting a new light on the question of what is the minimum requirement for a neuronal 
system to show criticality. To a large extent, the current literature has been focused on a form of homeostasis 
resulting from either a fine balance between excitation and inhibition, e.g., jlOllllj or some relatively complex 
dynamical processes at synaptic level, e.g., |15) . Our results show that a purely excitatory system can show 
the exact same behaviour such that on average each active neurone only activates one postsynaptic neurone. 
Here, this balanced state is achieved through a trade-off between the rates at which neurones become active 
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and quiescent. It should be noted that this formulation of the problem leads to interesting parallels with 
classical models of mathematical epidemiology which the authors intend to continue exploring. 

Spatial structure 

To make use of the analytic tractability of the mean field equation it was necessary to consider a fully 
connected network. While this is not true of the whole brain, it may be closer to the reality of the kind of in 
vitro systems typically considered in studies of neuronal avalanches. For example, Hellwig et al. j57] report 
up to 80% connection probability in local connectivity between pyramidal neurones in layers 2/3 of the rat 
visual cortex. Extending the work presented here to consider the effect of network topology on the system's 
dynamics and the resulting distribution of event sizes would be of particular interest from a developmental 
viewpoint. As networks mature, there is not only a switch to inhibition by a proportion of the neurones 
(the so-called GAB A switch), but also a subsequent pruning of synaptic connections [5S]. The level of 
pruning is high, with a 40% reduction in the number of synaptic connections between early childhood and 
adulthood ;58j. Thus, a developing network may be more readily approximated by a fully connected network 
than an adult neural network would be. 

The lack of a spatial embedding of our model is in contrast with many classical models of criticality, 
and also with physiological systems. Accordingly, our model cannot display another important marker of 
criticality, namely, the divergence of correlation lengths in space. A spatial embedding is not needed for our 
system to be critical and to exhibit a distribution of avalanche size similar to that observed in physiological 
neuronal avalanches. It therefore begs the question of the exact role of spatial embedding in the dynamics 
of neuronal avalanches. It may well be that, just like balanced activity in our model comes about from a 
trade-off between excitation and refractoriness rather than between excitation and inhibition, specific spatial 
embeddings may enable balanced activity without the need for plastic mechanisms. Kaiser and Hilgetag [59] 
showed that hierarchical modular networks can lead to limited sustained activity whereby the activity of 
neural populations in the network persists between the extremes of either quickly dying out or activating 
the whole network. Roxin and colleagues |60j observed self-sustained activity in excitable integrate-and-fire 
neurones in a small-world network, whose dynamics depends sensitively on the propagation velocity of the 
excitation. 
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Non-driven case 

Finally, in this paper, we have focused on the non-driven case h = 0. Whilst this constraint allowed the 
derivation of analytical results, it obviously contrasts with the reality of a physiological system unless one 
considers that any 'external' input operates at such a slower timescale that one could assume separation of 
time scales (an important assumption in the self-organised criticality framework). However, the fact that 
binning is required for identifying avalanches in physiological recordings suggests that this separation of time 
scales is unlikely. Whilst the introduction of a non-zero h in our model does not affect the results obtained 
using finite size expansion, it does effectively make it impossible for the system to operate at i?o = 1- A 
thorough investigation of the driven case (h > 0) will be the subject of the companion paper. 

Appendix - Van Kampen's system size expansion 

For generality we now assume that each neurone receives a constant external input and that the activation 
function can take forms other than the simple identity function. If Pn{t) is the probability that the number 
of neurones active at time t is n, then the master equation is given by 

^^=«(n + l)P„+i (t)- 
anPn (t) + 

f{^ + h) iN-n)P,,{t). 

Defining the step operator [e] (hereafter, in order to distinguish between operators and functions, we will 
use square brackets to denote operators and parentheses to denote functions, i.e., [^](/) means operator A 
acts on function /), 

M (.9 (")) =gin + l), 
[e-'] (ff(n))=g(n-l), 

we can now write 

= [s] (anP„(t)) - anP,,{t) + [e^i] (/ + h) {N - n)P„(t)) - f + h) [N - n)P„(t) 

= [{e - 1)] (anP„(t)) + [{s'^ 1)] (/ + h) (N - n)P„(i)). (18) 

The following analysis is based on introducing a continuous time-dependent density function, where 
the discrete state space becomes continuous. We do this by assuming that the fluctuations about the 
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microscopic value of n are of order . In other words, we expect that P„(i) wih have a maximum around 
the macroscopic value of n with a width of order N^. We thus take the ansatz that n(t) = Npi{t) +N^£^{t). 
Here fi{t) satisfies the macroscopic (mean field) equation and ^ is a stochastic perturbation scaled by 
. We can now write the distribution, Pn{t) as a function of ^, P„(t) ~ n(^,i). Doing so allows us to 
approximate the moments of the master equation using PDEs. 

Since [e] (n) = Nfi{t) + N^£_{t) + 1 = Nfi + N^£_ + N^^), we can say that [e] sends ^ ^ C + 
and we get 

[e] (P„(t)) = P„+i(t) ^ [e] (n(e,t)) ^m + N-i,t). 



Using Taylor's theorem we have 



[e] (/(0)-/(e + iV-5) = /(0 + 



1! dS, 2! 

Thus, wc can expand the step operator [e] (and similarly as 

\e] — I + iV^2 1 u 

e -^ =7— iV 2 1 

Now n = iV/i + A^^/^^ C " N^^n — fj.. Remembering we are interested in P„(t) for constant n wc then 
have 

dP„(t) _ dn(^, _ 9n dU _dU 1 rf^i 9n 
dt ~ di " at '^~dt^ 'dt ^ ^''dt'd^' 

Equating this to Equation 1181 we obtain 



N's-t- 

at dt 



' , 1 d iV-i 32 

\ 

9^ 2 de 



-N- 



.1 d iV-i 92 



2 9^2 

id 1 92 
,1 9 1 92 

-]Sf2 \ 

dc 2 9^2 



a[Nfi + N2q n 

/ (w (fi + N-in +h] (N- (Nfi + N^n ) n 



f(w(^i + N-H)+h) [l-(^l + N-^))n 
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For simplicity, we also now write f [w [ ii + N 2^] -\- h] as f. Thus, 



iV2 — 

dt dt 



' 1 d 1 92 

N2 ! . 

9e 2 de 
' 2de 



/ i-M-iv-^C n 



, 1 an amn a^l^^u 



2 oe 2 



^1 a 1 92 

-N2 1 

9^ 2 ae 

^1 a 1 32 

-N2 — \ 

a^ 2 ae 



-N2 



a 1 a^ 
a^ 2de 



(/n) 
(/A^n) 



an 



a^ 



2 ae 



N- 



ac 



,i5(n/) ^1 5(n/) 



9(n/o , ia2(n/) ^^^^(n/) iv-^a2(n/e) 



9^ 2 9^2 2 

We note we can expand / in terms of powers of N^'^ as follows 



ae 



a^ 



1 , w2£27V-l „ 

fiwii + h) + wN {w^l + h) + —^^^ / {w^i + h) + .... 



Writing / = /(^A* + h) and /' = f'{wfi + h) we find that 



9(/n) ^ an df 
a^ ■'as, as, 

^au i-'au 

= + N-i + nwf^ + (...) + .. . 

a^ifu) ^a^u 1 / ,-92n „ -:,aii\ , 



Similarly, 



g(/ng) ^ am , ,^^9/ 
9^ ■' as ^ ^'a^ 



= / 



9C 

5(110 

9^ 



dS 

dim 

9C 



new;/ )iv-i(...) + -. 
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Including these into our expansion we obtain 



dt dt ^ ae 2 

1 / -an 1 / -.an 

_ N2 f— + N-^ wif— + liwf 



ae V 

+ iv^A.(/f + iv-'(</'^ + n^/' 



u / -a^n 1 / -,a2n -,an 



+ N--- {■■■)■■ ■ 

1 an a(no and'^u ?an 
=^^"^aF + "^ + Tae"^'^ac 

- -a(no^/a2n -a^n 
+ 7v-^ (...).... 

We now group the terms in powers of A'^^: 

an ..la/ian ,ian/ - a a(ne) / -\ 
a^/aAi /_Ai/\ 

di"^ \ 2 2 2 j 

+ ^ (-</' + + ^nu-/' - nwf + iv-5 (. . .) + 
, a^n fafi f a(no 



ac^ \ 2 2 2 1 ac 

Choosing /i such that the macroscopic equation of the described process is satisfied (i.e. ^ = —a^i+{l — ^)f), 
we have the terms of order cancel out and we are left with an expansion in terms of iV^2 . The leading 
order is then 7V° (but relative to the macroscopic solution /i this is 0{N~^). nigher orders are now neglected 
as relative to the macroscopic ^ they are order (order of a single neurone). We are then left with the 
following linear Fokker-Planck equation 

an a(no / - a la^n 



dt ae (" + /-/'a-^)) + 2a^(-^ + (i-^)/) 
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It was shown [38j that the solution to this equation is gaussian so what we need to determine are the first 
and second moments. For ease we wih rewrite the above equation as 

a(no B{t) d^U 



an 

'dt 



A{t)- 



2 ' 

The derivative of the first moment with respect to time is then given by 



A 



B 
~2 



an 



Similarly, the derivative of the second moment with respect to time is given by 



di 



(19) 



dt 



R 9t 



,2 dm ,,,B f^^d'U 



A 



B 
"2 



9£, . 



= -2A{(,^)-B 
^~2A {i^)+B. 



(20) 



For both of these integrals many terms have vanished due to n(^, t) being a gaussian distribution and, in 
the limit, both the distribution and its derivative tend to zero quicker than any other factors involved. We 
can similarly calculate the variance of the fiuctuations as follows 



d{a^) _d{e) 



dt 



-2(0 



dt 



dt 

-2A{e)+B-2{(,){-A{0) 
-2^ [{e)-{i?^+B 



= -2A<T^ 



B. 



(21) 



Using Equations [19] and [21] with the correct values of A and B substituted in we can now compare the 
solutions of the following three equations with simulation results for the proposed model. We note that as 
long as we know the initial distribution of active neurones then for the initial conditions we have (Oo = 
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= ('^^)o ^ ^ ^^^^ (0 '^ill then remain zero for all time. 



|^ = -«M+(1-m)/, (22) 



d{0 
dt 



= -2 



a + / - - (0 , (23) 

(a + / - iz;/'(l - /i)) (^') + + (1 - m)/) ■ (24) 



at 

These in turn give the following equations for the mean, A, and variance, A^r, of the number of active 



neurones: 



A = N II + N ^ {^) = NjjL (assuming we know the initial number of active neurones), (25) 
A,^N{a^). (26) 
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Notes 

^As the power law is not a sufficient condition of criticality, one sliould not infer from this that the system is indeed critical, 
however, this step is commonly taken in published reports and that is worth mentioning here. 

^Strictly speaking the notion of critical slowing in neurones firing near firing threshold appeals to the same notion. 
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